NUMERICAL STUDY OF BREAKUP IN GENERALIZED 
KORTEWEG-DE VRIES AND KAWAHARA EQUATIONS 

B. DUBROVIN*, T. GRAVAt, AND C. KLEIN* 

Abstract. This article is concerned with a conjecture in [8j on the formation of dispersive 
shocks in a class of Hamiltonian dispersive regularizations of the quasilinear transport equation. The 
regularizations are characterized by two arbitrary functions of one variable, where the condition of 
integrability implies that one of these functions must not vanish. 

It is shown numerically for a large class of equations that the local behaviour of their solution 
near the point of gradient catastrophe for the transport equation is described locally by a special 
solution of a Painleve-type equation. This local description holds also for solutions to equations 
where blow up can occur in finite time. 

Furthermore, it is shown that a solution of the dispersive equations away from the point of 
gradient catastrophe is approximated by a solution of the transport equation with the same initial 
data, modulo terms of order where is the small dispersion parameter. Corrections up to order 
are obtained and tested numerically. 



Key words. Generalized Kortcweg-de Vries equations, Kawahara equations, dispersive shocks, 
multi-scales analysis 



AMS subject classifications. Primary, 65M70; Secondary, 65L05, 65M20 



1. Introduction. Many wave phenomena in dispersive media with negligible 
dissipation, in hydrodynamics, nonlinear optics, and plasma physics are described 
by nonlinear dispersive partial differential equations (PDE). These equations are also 
mathematically challenging since the solutions can have highly oscillatory regions and 
blowup even for smooth initial data. 

This article is concerned with a conjecture in [8] on the formation of dispersive 
shocks [16], [23], [26] in a class of Hamiltonian regularizations of the quasilinear trans- 
port equation 

ut + a{u)ux^0, a'{u)^0, u,xeM.. (1.1) 

In the present paper we will consider general Hamiltonian perturbations of (jl.ip up 
to fourth order in a small dispersion parameter < e <C 1. They can be written in 
the form of a conservation law 

ut + a{u)ux + e^dx {hi{u)uxx + &2(u)u^ + e \b3{u)uxxx + h{u)uxxUx + b'i{u)ul] 
[be{u)uxxxx + b7{u)uxxxUx + bsiu)ulx + bg{u)uxxul + bio{u)ul] } = 0, (1.2) 

where the coefficients 6i(u), . . . , bio{u) are smooth functions satisfying certain con- 
straints following from the existence of a Hamiltonian representation 

^ SH 
Ut + Ox . , , = 0, 
du(x) 

(see Corollary 12.21 below) . Here and below we use the notation 

Ox = TP- 

ox 
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This class of equations contains important equations as the Korteweg-de Vries (KdV) 
equation Ut + Guu^ + e^u^xx — and its generalizations, the Kawahara equation and 
the Camassa-Holm equation in an asymptotic sense (see [S]). 

Up to certain equivalencies the Hamiltonian regularizations of (|l.ip are charac- 
terized by two free functions c{u) and p{u): 



Ut + a{u)ux 



'dx 



ca'uxx + ^{ca'Yul 



+ €*dx y2pa' + -c^a" j Uxxxx + ■ ■ ■ = 0. 

(1-3) 

Two equations of the form (|1.3|) with the same invariants c(u) and p(u) commute, up 
to order ©(e^) 

where, for an arbitrary function a — a(u) 



ui + d{u)ux + e^dx 



cd'uxx + -^{ca'Yu^ 



+ e^dx 



2pa + -c^a" 
5 



0. 



In this paper the analysis of [8; up to order is extended to higher orders of e. 
It is shown that the only obstruction to the functions c{u) and p{u) by the condition 
of integrability is the condition that c{u) must not vanish. 

We then proceed to the study of the critical behaviour of solutions to (jl.2p . 
Namely, let (a^cj^cj^c) be a point of gradient catastrophe of a solution vP{x,t) to 
(jl.ip specified by an initial value vP{x,Q) = (j){x). This means that the solution is a 
smooth function of {x, t) for sufficiently small |a; — xd and t — tc<0. Moreover there 
exists the limit 



lim u {x, t) = Uc, 

but the derivatives u'^{x, t), (x, t) blow up at the point. The Universality Conjecture 
of [S] says that, up to shifts, Galilean transformations and rescalings, the behavior at 
the point of gradient catastrophe of a solution to (|1.2I) with the same e-independent 
initial data (j){x) essentially depends neither on the choice of the generic solution nor 
on the choice of the generic equation. Moreover, the generic solution near this point 
{xc,tc,Uc) is given by 



X — Xc — a(){t — tc) t — tc 



7 6 



4/7 



(1.4) 



where oq = a(uc) and the constants a, /3, 7 depend on the choice of the generic 
equation and the solution 

/126^ 



1 3 

■l22fc3(fe;)2 



1/7 



1/7 



(1.5) 



Here Oq — a'{uc), h\ = 61 (wc); it is assumed that h\ ^ Q. The constant k in these 
formula is inverse proportional to the "strength" of the breakup of the dispersionless 
solution vP{x, t) 

k — —6 lini =■ (1.6) 

where we assume that fc 7^ (another genericity hypothesis) and a'^k > 0. 

The function U = U{X;T), {X,T) e M^, is defined as the unique real smooth 
solution to the fourth order ODE 



X = TU - 



(1.7) 



which is the second member of the Painleve I hierarchy. We will call this equation 
PI2. 

The relevant solution is characterized by the asymptotic behavior 

U{X,T)^-{6X)i ^^^ + 0{X-i), X ^ ±00, (1.8) 

for each fixed T S M. The existence of a smooth solution of (|1.7p for all real X, T 
satisfying (|1.8p has been proved by Claeys and Vanlcssen 0. 

Observe that the principal term of the asymptotics (|1.4p depends only on the order 
regularization. In the present paper we will numerically analyze, in particular, the 
influence of the higher order correctionf0 on the local behaviour of solutions to (jl.2p 
near the point of catastrophe. 

First numerical tests of the PI2 asymptotic description of the critical point for the 
class of PDE in [8] have been presented in [T^ for the KdV and the Camassa-Holm 
equation. In [J] a rigourous proof of the asymptotic behaviour (|1.4I) has been obtained 
for the KdV equation. 

In this paper we generalize the numerical investigation of |15| to a larger class of 
equations which include the generalized KdV equation, the Kawahara equations with 
a dispersion of fifth order and the second equation in the KdV hierarchy. 

We comment on the formation of blow up and on the role of integrability in the 
formation of oscillatory regions. In particular we show the differences in the formation 
of dispersive shock waves between integrable and non-integrable cases. The KdV 
equation has been extensively studied numerically in |14| . 

Then we show numerically that the solution of the dispersive equation (|1.2p con- 
verges to the solution of the dispersionless equation (11.11) away from the point of 
gradient catastrophe at a rate of order e^. Finally we show that the solution of the 
dispersive equation (jl.2p is well approximated as a series in even power of e in terms 
of the solution of the dispersionless equation (II. ip up to order e'* by the so called 



^One should also take into account .8^ that the actual small parameter of the expansion 111.411 is 

(l2b;e2)l/^ 

In other words the asymptotic expansion II1.4I I makes sense only under the assumption 
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agreement with [2]. 



quasi-triviality transformation [5] away from the point of gradient catastrophe. Such 
an approximation has already been obtained for conservation laws with positive vis- 
cosity jl3J . Furthermore the existence of an expansion in even powers of e has already 
appeared and been proved in the context of large N expansions in Hermitian matrix 
models [3],[12]. 

The paper is organized as follows. In sect. 2 we briefly review the results of [5]. 
In sect. 3 we discuss higher order in e regularizations of (jl.ip and obstructions on the 
function c{u) by the condition of integrability. A numerical study of the applicability 
of the conjecture to generalized KdV equations is given in sect. 4. We also comment on 
the possibility of blowup. In sect. 5 the conjecture is tested numerically for equations 
with high order dispersion as the Kawahara equation. Differences in the formation 
of rapid oscillations in the solutions to integrable and non-integrable equations are 
studied. Details about the used numerical methods are given in the appendix. 

2. Hamiltonian PDEs and their invariants. In this paper we study scalar 
Hamiltonian PDEs of the order at most five. They are written in the form of a 
conservation law 

ut + dj:(p{u,eUx,e'^Ux^,e^Uxxx,e^Uxxxx) =0, (2.1) 

where 

Recall that the Euler-Lagrange derivative is defined by 

JIL^^_0^^ + d'^~ (2 3) 

Su{x) du ^ dux ^ duxx 

Here and in the sequel the integral of a differential polynomial is understood, in 
the spirit of formal calculus of variations^ as the equivalence class of the polynomial 
modulo the image of the operator of total x-derivative 

r. , dh dh 

Oxh = Ux^ + Uxx^ 1- (2.4) 

ou aux 

It is worthwhile to recall that a differential polynomial p{u; Ux, ■ ■ ■ , u^™-*) belongs to 
Imdx iff 

0, P ^ [ p{u;ux,...,u^"'^)dx. (2.5) 



5u{x 



The Poisson bracket of two local functional H, F associated with (|2.ip . (|2.2p . is a 
local functional of the form 

4 



Lemma 2.1. Equation (|2.1[) can be written in the Hamiltonian form (12. 2p iff the 
function if satisfies the following two constraints 



OUx 



dip 1 dip 
-A 



(2.7) 



Proof According to the classical Helmholtz criterion (see in [7]) the function (p{u, Ux, u^x, ■ ■ ■) 
can be locally represented as the variational derivative of some functional H — 
J h{u, Uxt ■ ■) dx iff it satisfies the following system of constraints 

5^ ^ y (!Ii±l)! 5^ ,^0,1,.... (2.8) 

m>0 

For the particular case under consideration the equations (|2.8p reduce to (|2.7p . □ 

Applying the Lemma to a PDE (|2.ip written in the form of the weak dispersion 
expansion one arrives at 

Corollary 2.2. The equation 

ut + a{u)ux + dx {ebo{u)ux + [biiu)uxx + b2{u)ul] + [b3{u)uxxx + hiu)uxxUx + b5{u)ul] 
[be{u)uxxxx + bj{u)uxxxUx + bsiu)ul^ + bg{u)uxxul + bio{u)u'^] } = (2.9) 

is Hamiltonian iff the coefficients bo, . . . , biQ satisfy 

bo = 0, 

^^3 = 0, 

b5 = Ib'^, 
b7 = 2b'e, 
bs = lb'„ 

bio = -^b'g. 

The Hamiltonian equations p.ip are considered modulo canonical transformations 
written in the form of a time-e shift 

u{x) ^ uix) = u{x) + e{u{x),K} + - {{u{x),K}, K} + . . . (2.10) 

generated by a Hamiltonian 

K ^ j k{u,eux,. . .)dx. (2.11) 

The transformations (|2.10p preserve the canonical form of the Poisson bracket (|2.6p . 
Two Hamiltonian equations are called equivalent if they are related by a canonical 
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transformation of the form (|2.10p . (|2.1ip . For example, the degree 3 terms in a 
Hamiltonian PDE of the form (|2.9I) can be ehminated by a transformation (|2.10p if 
a'{u) 7^ 0. Indeed, it suffices to choose the generating Hamihonian in the form 



K = e 



6a' (u) 

The following Lemma describes a normal form of Hamiltonians of order 4 (cf. [8]) 
with respect to transformations ()2.10|) . 

Lemma 2.3. Any Hamiltonian equation of the form (j2.9p with a'{u) ^ is 
equivalent to 



ut + a{u)ux + e^dx 



If 3 
hiu^x + -jj^Wx + [baUxxxx + ib'gUxxxUx + -^b'gul^ 



(2.12) 



0. 



The Hamiltonian PDEs \2.12\ and 



Ir 



biUxx + -xb'iUx + e bfiUxxxx + "^b'^UxxxUx + 7:b'eU^x 



(2.13) 



+bQUxxul + -b'gui 



0. 



are equivalent iff 



a = a, 6i = &i, ^6 = ^6 



(2.14) 



Proof We have already proved that the coefficients of degree 3 in e can be elimi- 
nated by a canonical transformation of the form (|2.10p . One can easily see that the 
coefficients a, bi and be are invariant with respect to these transformations. Two 
Hamiltonians of the form 



H 



£2 g4 



2 

XX 



12 



bgU^ 



dx 



and 



H 



f - —biu^ 



-beu 



2 

XX 



12 



bgU^ 



dx 



generating the flows (|2.12p and (|2.13p with the same coefficients a = a, bi = bi, 
be — be but with different 69 bg are related by a canonical transformation (j2.10p 
with 



K 



24 



■ u'tdx. 



□ 
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Thus the coefficients a, bi, bg, are invariants of the Hamiltonian PDE (I2.12[) . 

As it was discovered in [5] , any Hamiltonian PDE of the form (I2.12[) is integrable 
at the order approximation. More precisely, assuming a' ^ let us replace the 
invariants bi = bi{u) and bg = bQ{u) with 



61 k 
c = — , p 



a' ' ^ 2a' iU ' a'' 
Then the equation (|2.12p is equivalent to the PDE 



ut + a{u)ux + e^d; 
with the Hamiltonian 



ca'uxx + ^{ca'Yul 



2pa' + gC^a" ) Ux 



r - rill 2 I 4 



( till I 3 2 r(i)\ 2 



(2.15) 



= 
(2.16) 



(2.17) 



dx 



where, as above, 



nu) = a{u). 



The approximate integrability means that, fixing the functional parameters c = c(w), 
p — p{u) one obtains a family of Hamiltonians satisfying 



(2.18) 



for an arbitrary pair of smooth functions / — J{u), g — g{u). In particular choosing 
f[u) — hu^ one obtains the Hamiltonian 



H 



2 2 I i I \ 2 



dx 



6 2 

of a general order 4 dispersive regularization of the Hopf equation 
1 



(2.19) 



ut+uux+e dx 



cuxx + T:c'ui 



introduced in [5] 4 

More generally, we call a perturbation 



+e^dx \2puxxxx + ^p'uxxxUx + "ipulx + 2p"uxxU^x] = 

(2.20) 



of the Hopf Hamiltonian 



Hf) = — dx 



6 



^In the present paper we use a different normalization c(u) 1— > 12c(u). 
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N -integrable if, for any smooth function / — f{u) there exists a perturbed Hamiho- 



fe>i 



such that for f = \ the Hamiltonian Hf coincides with H and, moreover, for any 
pair of functions /, g the Hamihonians Hf, Hg satisfy 



{Hf,Hg)^0{e^+'). 

For example, the perturbed Hamihonian (I2.19P is 5-integrable. The commuting 
Hamiltonians have the form (I2.17p . In the next section we wiU discuss the problem 
of constructing higher integrable perturbations of (|2.19p . 

3. On obstacles to integrability. We will now study the possibility to extend 
the commuting Hamiltonians (|2.17p to the next order of the perturbative expansion. 

Theorem 3.1. 1) Any order 6 perturbation of the cubic Hamiltonian Hq = 
J ^ dx can be represented in the form 



H 



6 



-c{u)ul + e*p{u)ul^ - {a{u)ul^^ + 



dx 



(3.1) 



Such a perturbation is 7-integrable for arbitrary functional parameters c — c{u), p 
p{u), a = a{u), P = I3{u). 

2) The perturbation (|3.ip can be extended to a 9-integrable one iff c{u) ^ and 



a 



28 



67p c' + 33 cp' + 12 cc'^- 9 c" 



(3.2) 



Proof A general order 6 perturbation of the cubic Hamiltonian Hq must have the 
form 



H = 



The last two terms can be eliminated by a canonical transformation 

H^^ H ~ t{H,F} + ... 

with 



1 



For an arbitrary function / = f{u) the density of a Hamiltonian 



H 



f 



hf dx 



commuting with p.ip modulo O (e^) must have the form 

[pf" + ^c'f^'^) ul^ - \ (p'f^^) + ^cc"/(^) +Pf^'^ 

(3.3) 

[af{u)ul^^ + I3f{u)ul^ + 7/(m)m|^u^ + 5f{u)ul] 



hf^f-^-cf"ul + e' 



with some smooth functions af = af{u), /3f = (3f{u), jf — 7/(u), 6f — 5f{u) 
depending on /. From the commutativity 



{H,Hf} = 0{e') 
one uniquely determines these coefficients 



af = a 



'3 253pc' + 169cp' cc'^ 5 2 A aa) /^29 31 



'^■^ '^■^ [^2 168 35 56 ^ ^21 70 / 

7/ = - y"' + - c'c'" - 3cc'c") + c'p' - ^pc" - cp") Z^^) 

pc' + |(cc'^ + Uc^c") + ^cp') /(^) ~ ^ (23cp + gc^c') /f^) - ^c'fi^) 



37 



/ 1 , , 10CCV" + 7CC^CW+CV5) 2 1 (4A 

1 „ 1 , „ 3 „2 3 „, cc'c'" 1 ,„ c2c(4)\ , . 

— a" + -p'c" H cc" H pc'" H 1 cp'" H /^^^ 

15 5^ 40 10^ 10 15 ^ 15 y 

+ + 15^^ + 3^^ + 40 + 10'^^ j ' 

+ (j_a + -pc' + — + -cp' + -c^c"] /(^) +(-cp+ -c^A /(«) + — 
15 6^ 16 10 ^ 40 J ■' V20 80 J ■' 240'' 

Thus the resulting Hamiltonian Hf satisfies 

It is not difficult to also verify the commutativity 

{Hf,Hg}^0{e^) 

for an arbitrary pair of functions / — f(u), g = g{u). 

Let us now analyze the possibility of extension to a commutative family of order 
8. We add to p.ip terms of the form 



4 

XX 



+ ^5W^U^2.j. + AqUxxUxxx + A7Uxxxx] 

and to (|3.3p a similar expression 

+^5MlwLa; + Bqu^^uI^^ + Brul^^^] dx. 



4 

XX 



Here Ai, A^, Bi, Bj are some functions of u. The goal is to meet the 
condition 



{H,Hf}^0{e'). (3.4) 

The order 8 terms in the bracket p.4p are represented by a differential polynomial of 
degree 9. From the vanishing of the coefficient of u'^^'>Ux it follows that 

Next, from the vanishing of the coefficient of u^^^^xx we get (|3.2p . 

Further calculations allow one to determine Bq from the coefficient of u^^^u^xUx, 
Bq from the coefficient of i?4 from the coefficient of u^^^ u^^Ux , B3 from the 

coefficient of u'^^^UxxU^, -B2 from the coefficient of u^^^u^, and, finally, Bi from the 
coefficient of Uxxul- All these coefficients are represented by linear differential oper- 
ators of order at most 12 acting on the arbitrary function / = f{u). The coefficients 
of these operators depend linearly on , . . . , and their u-derivatives and also on 
c(u) and p(u) and their derivatives. The explicit formulae are rather long; they will 
not be given here. As above one can verify validity of the identity 

for any pair of functions /(u), g{u). □ 



Corollary 3.2. Let p{u) be an arbitrary non-vanishing function. Then the 
Hamiltonian 



H = 



y + e p{u)u 



dx 



(3.5) 



cannot be included into a 9-integrable family. 

4. Quastriviality transformations and perturbative solutions. In this 

section we will develop a perturbative technique for constructing monotone solutions 
to the equations of the form (I2.20p for sufficiently small time t. This technique is 
based on the so-called quasitriviality transformation [S] expressing solutions to the 
perturbed equation (j2.20p in terms of solutions to the unperturbed equation. 
To explain the basic idea let us consider the equation 



ut+uux + e'^dx ( cuxx + \cul ) H =ut + dx 



H = 



2" + 



dx. 



5H 

' Su{x) 



= 



(4.1) 



The quasitriviality transformation for this equation 



V ^ u ■ 
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+ c'Vxx + ^c"vl 



(4.2) 



is generated by the Hamiltonian 



K 



(4.3) 



u^v + e{v{x),K} + -^{{v{x), K}, K} + 



Substituting into eq. ()4.ip one obtains a function u{x^ t; e) satisfying (|4.ip up to terms 
of order e^. Indeed, one can easily derive the foUowing expression for the discrepancy 



uux + e^dx { cuxx + ^c' 



Ut 



'^'^^xx Vxxx , "^^^xx ^xxxx , 57 Vxx'^xxx ^^xx^xxxxx ^^'^xxx'^x 



2vx 



4 Vrr. 



4Vx^ 



+ 



-'XXXXXX 



+cc' 



+CC 



Avx^ 



2,Vx 



+ 



1 3 XXX 



4w^ ^ 2~ 



2W2; 

2W2: 



+ 3 VxxV. 



XX ^xxx 



+ CC 



a: ^xxx 



■ SVxVxx 
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2 Vx 



'^XXX \ , I n ( 77; 

+ C C 



a; ^xxx 



lOvxVxx"^ + 9c'c"'t;3„^^ + c'c^^^l 



+^c"\lvxx + -/C"'vl 



Note that the same quasitriviality transformation works for solutions v = v(x,t) to 
the nonlinear transport equation 

vt + a{v)vx = 

transforming it to solutions, modulo 0(6^), to the perturbed equation (II. ip . 

Denote by 4'{x) = v{x, 0) the initial data for the Hopf equation. The initial value 
of solution u(x,t;e) given by the formula (|4.2I) differs from (j){x): 



u{x, 0; e) = (f) + 



+ CI(e'), (4.4) 



In order to solve the Cauchy problem for (14. ip with the same initial data one can use 
the following trick. Let us consider the solution v — v{x,t]e) to the Hopf equation 
with the e-dependent initial data 



v{x, 0; e) ^ (f) — 



Vxx \ I / / 
— ^ + C 0: 



XX 2 



Such a solution can be represented in the form 

v{x, t; e) = v{x, t) + e^w{x, t) + ©(e") 
where the function w{x,t) has to be determined from the equation 



^'\v) 2$'(w) 



(4.5) 



(4.6) 



(4.7) 
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Here is the function inverse to 4'{x). Applying the quasitriviahty transformation 
to the solution v(x, t; e) one obtains a function 



u{x, t;e) = v + e^w + 



' + c'v^^ + -d'vi 



(4.8) 



2\ vl ] ^^2 

satisfying equation (|4.ip modulo terms of order e^* with the initial data 

u(a;,0;e) = 0(a;)+O(e4). (4.9) 

5. Generalized KdV equations. In this section we will first study the role 
of the function a(u) in (II. ip on the validity of the conjecture. This is done for the 
generalized KdV equations having the form 

Mt + a(-u)-U:E + e^Wxxa; = 0. (5.1) 

We will assume that a(u) is monotonic in an open neighborhood of each critical point. 
The functional parameters c{u) and p(u) in (I2.17|) are given by 

The basic idea of the PI2 approach to the breakup behavior is that the equation 
behaves in this case approximately as the KdV equation. We will test this assumption 
first for a(u) of the form a(u) = 6u", n S N. 

5.1. Breakup. To begin we will study the solutions to generalized KdV equa- 
tions close to the breakup of the corresponding dispersionless equation. A generic 
critical point {xc^tc^Uc) is given by 

a{uc)tc + $(Uc) = Xc, 
a'(u,)t, + $'(u,) = 0, (5.3) 
a"{u,)t, + $"(u,) = 0, 

where $(u) is the inverse of the initial data (j){x) (which might consist of several 
branches). We will always study the initial data ^(a;) = sech^a;, which imply $(m) = 
ln((l ± ^/l~^)/^/u). For the critical values we obtain 

2n 



2n + 1 ' 
( l + 2n)" +i/^ 

V2n+1 , / V2n + 1 + 1 
'In 



2n V \/2^ 

k = -i(a"'(u,)i, + <i>"'(^.,)) - ^^"^ + 



6 96n2 

We study first the difference between the numerical solution to the generalized KdV 
equation and the solution to the dispersionless equation, a generalized Hopf equation, 
on the whole computational domain. For values of e = 10^^, 10^^'^^, . . . , 10"'^ we find 
that this difference scales for n = 1 (KdV) roughly as e" with a — 0.299 (correlation 
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coefficient r = 0.99997 in linear regression, standard deviation — 0.0018), for 
n = 3 we have a = 0.317 (r = 0.9998, = 0.0046), for n = 4 we liave a = 0.324 
(r = 0.9998, = 0.005), and for n = 5 we have a = 0.325 (r = 0.9998, = 0.0053). 
The predicted value is 2/7 = 0.2857. It can be seen that the above values are all 
higher, and that the scaling for the generalized KdV equations is close to e^^^. Thus 
the decrease is at least of the predicted order. It is not surpising that higher values 
for the exponent are found since we consider considerably large values of e for which 
the contributions of higher order in the difference still play a considerable role. 

As discussed in the previous section it is conjectured that the behavior of the 
solutions to the generalized KdV equation in the vicinity of the critical point is given 
in terms of the special solution to the PI2 equation. Expanding a{u) for u ^ Uc as in 
[8], one finds the behavior shown in Fig. 15.11 for different values of n. It can be seen 
that the asymptotic description is much better for KdV due the fact that the PI2 
transcendent gives an exact solution to KdV. For other values of n this transcendent 
gives the conjectured description in the vicinity of the critical point. The quality of 



n=1 



n=3 





0.75 0.8 0.85 0.9 

X 

n=5 




0.55 0.6 



Fig. 5.1. The blue line is the solution of the generalized KdV equation ut + (}u"'Ux + e^Uxxx = 
for different values of n for the initial data cl>(x) = 1/cosh'^a; and e = 10~^ at the time tc and near 
the point of gradient catastrophe Xc of the Hopf solution ( center of the figure). The green line is the 
multiscale approximation in terms of the PI2 solution. 

the asymptotic description shows the expected scaling for smaller values of e as can 
be seen on the left in Fig. 15.21 

The quality of this PI2 approximation is not limited to functions a{u) in ()5.ip 
which are polynomial in u. If we consider the case a(u) = 6sinhu, we obtain the right 
figure in Fig. 15.21 It can be seen that the PI2 asymptotics gives the same excellent 

13 



0.634 0.636 0.638 0.64 0.642 0.644 0.646 0.64B 0.65 0.652 



1.3 1.32 1.34 1.36 1.38 1.4 1.42 1.44 1.46 1.48 



Fig. 5.2. The blue line is the solution of ut + a(u)ux + e^Uxxx = with a{u) = 6u^ on the left 
and a{u) = 6sinhn on the right, for the initial data (f>{x) = 1/ cosh^ x and e = 10~* at the time tc 
near the point of gradient catastrophe Xc of the Hopf solution ( center of each figure). The green line 
is the multiscale approximation in terms of the PI2 solution. 



description as for KdV. 

5.2. Oscillatory regimes and blowup. It is known that solutions to initial 
value problems with sufficiently smooth initial data for the generalized KdV equation 
with n < 4 are globally regular in time. This is not the case for for n > 4 where 
blowup can occur at finite time with n = 4 being the critical case. For this case a 
theorem by Martel and Merle [24] states that solutions on the real line, with negative 
energy, blow up in finite or infinite time. For the general case n > 4 and periodic 
settings considered here, the question is still open. Since the energy has the form 



E 



7/^ 

2 " (n + l)(n + 2) 



dx 



it will be always negative for sufficiently small e and positive u. 

Here we address numerically the question whether the formation of dispersive 
shocks, i.e., of a region of rapid modulated oscillations, precedes a potential blowup. 
We expect that the breakup of the solution to the dispersionless equation is regularized 
by the dispersion in the form of oscillations which then develop into blowup if the 
latter exists. This is exactly what we see in the following. Notice that the breakup 
time is given by the dispersionless equation and is thus independent of e. We first 
study the case n = 4. For e = 1, the energy is positive and no indication of blowup is 
observed. For e = 0.1 we obtain for the initial data 0(a;) = sech^a; the left figure in 
Fig. 15.31 For smaller e (e = 0.01) the behavior is similar, but there are as expected 
more oscillations, and the size of the oscillations reaches higher values earlier as can 
be seen in Fig. 15.31 For obvious reasons it is numerically difficult to decide whether 
the dispersive shock will lead to a blowup. In practice we run out of resolution before 
the code breaks down because of a blowup. This is due to oscillations in Fourier space 
as can be seen in Fig. 15.41 Though there is in principal enough resolution to approach 

t), the oscillations of the Fourier coefficients make an accurate approximation via 
a Fourier transform impossible. 

The reason for this behavior is as discussed in f27| that singularities of the form 
(z—Zj)^^ in the complex plane lead asymptotically to Fourier coefficients with modulus 
of the form Cfc~('^J+^) exp(— (5/c), (5 > 0. If there are several such singularities, there 
will be oscillations in the Fourier coefficients. In the present case there are at least 
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Fig. 5.3. Solution of the generalized KdV equation with n = 4 for the initial data 4>{x) = 
1/ cosh^ X and e = 10~^ at the time t = 0.3180 ^ tc on the left, and for e = 10~^ at the time 
t = 0.2235 > tc on the right. 




0.35 



0.05 0.1 0.15 0.2 0.25 

t 

Fig. 5.4. Logarithm with base 10 of the 
modulus of the Fourier coefficents of the func- PlG. 5.5. L°° norm of the solutions in 

tion shown in Fig. 15.31 on the right. pig^ \5.S[ 



two such singularities, the breakup which is strictly speaking only singular for e = 0, 
but has an effect already for finite e, and the blowup, which leads to the behavior 
seen in (|5.4p . In Fig. 15.51 we give the L°°-norm of the solutions in Fig. 15.31 It cannot 
be decided on the base of these numerical data whether there is finite time blowup in 
this case. If it exists it is clearly preceded by a dispersive shock. 

In the supercritical case n = 5 we obtain a similar picture. In Fig. I5.6l we see the 
solution in the case e = 0.1. Again it appears as if the rightmost peak evolves into a 
singularity. For smaller e (e = 0.01) there are again more oscillations, which stresses 
the importance of dispersive regularization before a potential blowup. Studying the 
i°°-norm of the solutions in Fig. 15. 6[ we can see that the case e = 0.1 indeed seems 
to approach an L°° blowup in finite time. Because of resolution problems we could 
not reach a similar point for e = 0.01. 



6. Kav^rahara equations. The Kawahara equations |20) which appear in general 
dispersive media where the effects of the third order derivative is weak as in certain 
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Fig. 5.6. Solution of the generalized KdV equation with n = 5 for the initial data (j){x) = 
1/ cosh^ X and e = 10~^ at the time t = 0.2362 3> tc on the left and for e = 10~^ at the time 
t = 0.2235 > tc on the right. 




0.12 0.14 0.16 



Fig. 5.7. L°° norm of the solutions in Fig. 15.61 



hydrodynamic or magneto-hydrodynamic settings can be written in the form 

ut + ^dxfiu, e Ux, e^u^x) + P e'^Uxxxxx = 0. (6.1) 
Here we will mainly study the case 

f{u,eux,e^uxx) = + '2ae^Uxx (6.2) 



with a = 1 and 13 — ±1. The global well posedness of solutions of (|6.ip in a suitable 
Sobolev space has been proved in |25j . 

The functional parameters c{u), p{u) in (|2.17p are constants 

c(m) = ^a, p{u) ^ -^(3. (6.3) 

At the critical point we obtain for /3 = — 1 that the breakup behavior is well 
described by PI2 in lowest order as can be seen in Fig. 16.11 
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Fig. 6.1. The blue line is the solution of the Kawahara equation (f) = —I) for the initial data 
(t){x) = 1/ cosh^ X and t = 10~^ on the left figure and t = 10~'^ on the right figure. The plot is taken 
at the time tc near the point of gradient catastrophe Xc of the Hopf solution (center of the figure). 
Here Xc ~ 1.524, tc ~ 0.216. The red line is the corresponding Hopf solution, the green line the 
multiscale approximation in terms of the PI2 solution. 

It can be seen that the PI2 sohition gives close to the breakup point a much better 
description of the Kawahara solution than the corresponding Hopf solution. The 
oscillation closest to the breakup point is too far away from the latter to be correctly 
reproduced, but the PI2 solution catches qualitatively the oscillatory behavior of the 
Kawahara solution near the critical point. With smaller e, the agreement gets as 
expected better, see the right figure in Fig. 16.11 




1.44 1.46 1.48 1.5 1.52 1.54 1.56 1.58 1.6 1.62 1.44 1.46 1.48 1.5 1.52 1.54 1.56 1.58 1.6 1.62 



Fig. 6.2. The blue line is the solution of the Kawahara equation (13 = 1) for the initial data 
(f>{x) = 1/ cosh^ X and e = 10~^ on the left figure and t = 10~^ on the right figure. The plot is taken 
at the time tc near the point of gradient catastrophe Xc of the Hopf solution (center of the figure). 
Here Xc ~ 1.524, tc ~ 0.216. The red line is the corresponding Hopf solution, the green line is the 
multiscale approximation in terms of the PI2 solution. 

For (3 — 1, the breakup behavior of solutions to the Kawahara changes as can be 
seen from Fig. 16.21 In this case the oscillations in the Kawahara solution appear on 
the other side of the critical point and around it with small amplitude. This behaviour 
cannot be captured by the PI2 solution, but it is a higher order effect. Close to the 
critical point, the multiscale solution gives as before a much better description of the 
Kawahara solution than the Hopf solution. 
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For smaller values of e, both asymptotic solutions become more satisfactory as 
can be seen from the right figure in Fig. 16.21 It is interesting to notice that with 

E=1 0'^ 

1.5 
1 

= 0.5 


-0.5 

1.6 1.8 2 2.2 2.4 2.6 

X 

E=10"^ 



liiiiiiii^^ 



1.6 1.8 2 2.2 2.4 2.6 

X 

Fig. 6.3. Oscillatory zone of the solutions to the Kawahara equation ((3 = 1) for the initial 
data </)(a:) = 1/ cosh^ x and two values of e at time t = 0.25 > tc- 

decreasing e, the oscillations become smaller in amplitude in this case, but appear 
closer to the critical point. It can also be seen that the solution has the tendency to 
form one oscillation on the other side of the critical point close to the corresponding 
PI2 oscillation. Tracing the solution for larger times, it can be recognized that this 
will be the only oscillation to the left of the critical point, whereas a zone of high- 
frequent oscillations which appears to be essentially unbounded (see [H]) develops to 
the right, see Fig. 16.31 The oscillations appear to be as in the KdV case more and more 
confined to a zone similar to the Whitham zone, though no asymptotic description of 
the oscillations exists since the equation is not integrable. 

It seems also that this one oscillation to the left is really due to the third order 
derivative in the Kawahara equation as can be seen in Fig. 16.41 where one has to the 
left the Kawahara solution from Fig. 16.21 and to the right the analogous solution for 
a = 0, i.e., Kawahara without third order derivative. The oscillations to the right of 
the critical point being due to the fifth order derivative are present in both cases and 
have only slightly different form. 

6.1. PDE with nonlinear dispersion. To show that the breakup behavior 
discussed in the previous sections is typical, we will now consider equations of the 
form ()2.16p with nonlinear dispersion, i.e., with functions c{u) and p{u) not constant. 
The Camassa-Holm equation (CH) falls in this class if the nonlocal term is expanded in 
a von Neumann series, see P|, for functions c(u) ~ u and p{u) ~ u. The applicability 
of the PI2 asymptotics to CH was studied numerically in [15]. For simplicity we 
restrict our analysis to the case c{u) and p{u) both proportional to with the Hopf 
equation as the dispersionless equation and initial data of the form (t){x) — sech^x. 
More complicated functions c and p can be considered, but the results are qualitatively 
the same. 

In Fig. l6.5l the behavior at the critical time can be seen for c(u) = and p{u) = 0. 
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Fig. 6.4. Solutions of the Kawahara equation \6.1\ with f) = 1 a = 1 to the left and a = to the 
right. The solution is given for the initial data <f>(x) = 1/ cosh^ x and e = 10"'^ at the time tc near 
the point of gradient catastrophe Xc of the Hopf solution. 



The situation is obviously as in the KdV case. 




Fig. 6.5. Solution to the equation i2.16V 
for a{u) = u, c{u) = andp(u) = and ini- 
tial data i}>{x) = sech?x at the critical time, 
and the corresponding multiscale solution in 
terms of the PI2 transcendent. 



Fig. 6.6. Solution to the equation i2.16V 
for a{u) = u, c{u) = —p{u) = and ini- 
tial data <j){x) = sech^x at the critical time, 
and the corresponding multiscale solution in 
terms of the PI2 transcendent. 



As for the Kawahara equation the relative sign between the third and the fifth 
derivative is important for the form of the oscillations. The situation with the opposite 
sign of c and p can be seen in Fig. 16.61 It is qualitatively the same as in the KdV case. 
New features appear as in the case of the Kawahara equation for the same sign in 
front of the third and fifth derivative. As can be seen in Fig. 16. 7[ oscillations of small 
amplitude appear as in the Kawahara equation on the other side of the critical point. 
Thus non constant functions c{u) and p{u) as expected do not change the picture 
from the case of constant functions as long as they do not vanish at the critical point. 

6.2. Quasi-trivial transformation. In this subsection we study numerically 
the validity of the expansion given in sect. 4. For times t ^ tc the behavior of 
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Fig. 6.7. Solution to the equation f2. 1 6\) for a{u) = u, c{u) = p(u) = and initial data 
(f>{x) = sech^x at the critical point, and the corresponding multiscale solution in terms of the PI2 
transcendent. 

the solution of (|4.1[) should be described to order by the solution of the Hopf 
equation ut + uu^ — with the same initial data. In fact we find that the difference 
between the Hopf solution and the solution to Kawahara equation (|6.ip with /3 = 1, 
a = 1 for the initial data (t>{x) = sech^a; at t = ic/2 scales as with 7 = 1.94 
(values of e = 10^^, 10^^'^^^, . . . , 10~'^, correlation coefficient r = 0.9997 in linear 
regression, standard deviation CTq. = 0.027). For the same setting in the interval 
X S [0.8, 2], the difference between the quasitriviality solution as described in sect. 4 
and the Kawahara solution scales as with 7 = 3.77 (correlation coefficient r = 0.999 
in linear regression, standard deviation Ua = 0.088). This confirms the theoretical 
expectations. In Fig. 16.81 the difference between the Kawahara and the Hopf solution 
and the quasitriviality transform in order can be seen for this case. A similar 
scaling is observed for a = —1 and a = (KdV). 

The difference between the solution of generalized KdV equation ut + u^u^ + 
^'^Uxxx — and the solution of the corresponding conservation law scales, dX t — tc/2 
as with 7 = 1.9890 (values of e = 10~^, 10~^-^^^, . . . , 10"'^, correlation coefficient 
r = 0.99998 in hnear regression, standard deviation CTq = 0.0068). 

6.3. Second equation in the KdV hierarchy. In this subsection we study the 
formation of dispersive shock waves for a family of equations that includes integrable 
and non- integrable PDEs. Interestingly the family of equations 

Ut + iQu^Ux + 10ae^(uu 

XXX ~t~" '^'^x'^X'X^ ~t~ C Uxxxxx — 

(6.4) 

having the invariants 

cm = —a, p{u) = 



is completely integrable for a — ±1 and coincides with the second equation in the KdV 
hierarchy (KdVII). Varying this factor one can study the transition to the Kawahara 
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Fig. 6.8. Solution to the Kawahara equation with a = 1 for the initial data 4>{x) = 1/ cosh'^ a; 
at the time t = tc/2 for two values of e; in blue the difference between the Kawahara and the 
corresponding Hopf solution, in green the order term of the quasitriviality transformation. 



equation. As expected KdVII shows similar oscillations as KdV [Tl], as can be seen 
Fig.lOl 

For larger values of a one can recognize in Fig. l6.9l also a formation of oscillations 
on the other side of the inflection point as in the Kawahara equation. These effects 
become smaller for larger values of a (a > 1.2), but it shows that the phenomenon of 
integrability with the appearance of KdV-type oscillations is rather subtle. Thus it 
seems that the decisive factor for the appearence and the size of these oscillations is 
the relative sign and size of the factors in front of the third and the fifth derivative in 
the equation. Notice that equation (16.41) is for smaller a closer to the non-integrable 
(in higher orders in e) equation 

(see Section [3] above) . It would be interesting to elaborate this observation in order 
to develop numerical tests of (approximate) integrability based on the study of the 
phase transition from regular to oscillatory behaviour. 
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Appendix A. Numerical Methods. In this appendix we will briefly review 
the used methods in the numerical study of the PDE in the small dispersion limit and 
of the PI2 solution and give references in which details can be found. 

Since critical phenomena are generally believed to be independent on specific 
boundary conditions, we restrict our analysis to essentially periodic functions. Typi- 
cally we consider Schwarzian functions on a domain on which the functions are at the 
boundaries smaller than machine precision (10~^^ in double precision). Such func- 
tions can be periodically continued and are smooth with numerical precision. This 
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Fig. 6.9. Oscillatory part of the solution to the KdVII equation for the initial data (t>{x) = 
1/ cosh^ X and e = 10~^ at a time t = 0.04 > tc = 0.029 for several values of a. 



allows a Fourier discretization of the spatial variables and an approximation of the 
solutions via truncated Fourier series. The use of Fourier spectral methods is espe- 
cially efhcient for the studied dispersive PDF because of the excellent approximation 
of smooth functions and the only minimal introduction of numerical dissipation. The 
latter is especially important if one is interested in the study of dispersive effects. 

After discretization of the spatial coordinates, the PDF is equivalent to a typically 
large system of ordinary differential equations (ODF) in the time variable. Because 
of the high order of the spatial derivatives and because of the strong gradients we 
want to study, these systems will be typically stiff. If the stiff part is linear as is the 
case for the generalized KdV equations and for the Kawahara equations, the system 
of ODF has the form 

Lv + N[v] = 0, 

where v is the discrete Fourier transform of the solution, where L is the stiff linear 
operator, and where the nonlinear term N[v] contains only derivatives of lower order. 
For such systems, efficient integration schemes exist. We use a fourth order exponen- 
tial time differencing scheme [6], see [21j for a comparison of fourth order schemes 
for KdV. The numerical accuracy is controlled by sufficient spatial resolution, i.e., 
Fourier coefficients decreasing to at least 10~*, and by numerically checking energy 
conservation. Since all equations studied here are Hamiltonian, energy is a conserved 
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quantity. Due to unavoidable numerical errors, it will be weakly time dependent in 
numerical time integrations. As discussed in |21j . conservation of the numerically 
computed energy typically overestimates the accuracy of a solution by two orders of 
magnitude. We always compute with an error in energy conservation smaller than 
10^^ which implies that the error is well below plotting accuracy. 

The situation is different for the equations with nonlinear dispersion in sect. 6.2. 
For these PDE we use an implicit fourth order Runge-Kutta method (Hammer and 
Hollingsworth method). These equations are numerically much more demanding. 
Therefore we compute with lower spatial resolution and an energy conservation of the 
order of lO""*. 

The special solution of the PI2 equation is generated with the code bvp4 dis- 
tributed with Matlab. For details see [15]. The Hopf solution is obtained from the 
implicit form u{x,t) = (j^i^)-: ^ = t4'{0 + ? with a fixed point iteration to machine 
precision. The derivatives of the Hopf solution are obtained by evaluating the the 
analytic expressions following from the characteristic method. 
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